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Abstract 

Due to the teed back of the user community, 
three major features have been added to the NASA 
Lewis ice accretion code LEWICE. These features 
include: first, further improvements to the numerics of 
the code so that more time steps can be run and so 
that the code is more stable; second, inclusion and 
refinement of the roughness prediction model 
described in an earlier paper; third, inclusion of multi- 
element trajectory and ice accretion capabilities to 
LEWICE. This paper will describe each of these 
advancements in full and make comparisons with the 
experimental data available. Further refinement of 
these features and inclusion of additional features will 
be performed as more feedback is received. 


Nomenclature 

b = bead height/roughness (m) 
F = wetness fraction 
V = velocity (m/s) 
s = surface distance (m) 
r = water density 
s = surface tension (kg/m/s) 


Introduction 

Three major features have been added to the 
NASA Levins ice accretion code LEWICE . 1 The 
numerics have been improved, the roughness predic- 
tion model has been revised, and a multi-element 
capability has been added. 

The numerics of the code have been improved in 
several ways to produce more accurate ice shapes. 
The convergence characteristics of the code have 
been improved by implementing an adaptive grid 
technique, a new ice growth algorithm and a new 
variable time stepping scheme. Improvements to the 
transition model and transition heat transfer calcula- 
tions have been made to produce more realistic 
results. New additions include a “pseudo" surface 
which produces more realistic heat transfer for large 
glaze ice shapes and a mass addition routine which 
allows ice growth in arbitrary directions. 

An adaptive grid scheme has been implemented, 
which allows more optimal tailoring of the individual 
surface models for each phase of the ice growth pro- 
cess yielding smoother more accurate ice shapes, 
better convergence characteristics and quicker run 
times. A highly refined “baseline” model is used to 
represent the geometry at each time step. This model 
is updated after each time step. 

Adaptive griding techniques are used to generate 
optimal surface models from the baseline model for 
the heat transfer, the collection efficiency, mass bal- 



ance, energy balance, and the mass addition 
phases. A typical baseline model may contain 4000 
points. Atypical flow panel model for accurate trajec- 
tory calculation may have constant leading edge 
spacing and require 100 points. Flow panel models 
used for the generation of velocities for the heat 
transfer coefficient calculation, which can be either of 
the “pseudo” surface type (which produces more 
realistic heat transfer for glaze shapes) or keyed to 
radius of curvature, may typically require 150 panels. 
A typical surface model for the energy balance, the 
mass balance and mass addition algorithms may 
contain several thousand points. 

A new ice growth scheme has been imple- 
mented. This scheme employs a separate time step- 
ping procedure on a highly refined surface model. 
The method features local conservation of mass, 
accurate resolution of complex ice shapes, and erad- 
ication of the troublesome problem of iced lobes 
growing into each other. 

The time stepping procedure has been auto- 
mated. The user specifies the maximum ice thick- 
ness to be added at each time step. The time step is 
varied at each time step to match this maximum 
amount of ice thickness. This ice growth scheme 
yields better convergence characteristics by giving 
the user more precise control over the maximum 
geometric change and hence aerodynamic change 
between consecutive time steps. 

Additional improvements include improvements 
to the transition model, the transition heat transfer 
calculation, incorporation of a more realistic “pseudo" 
surface and additions to the ice growth module. The 
transition model has been modified to produce more 
reafistic transition locations for cases with large or 
multiple stagnation points. A more realistic treatment 
of the laminar and turbulent heat transfer coefficient 
in the transition region has been implemented. The 
new ice growth model allows ice growth in arbitrary 
directions to accommodate current and future ice 
growth models. An optional “pseudo" surface method 
has been installed which more accurately models 
flow and hence heat transfer for ice shapes with 
large stagnation zones. 

The roughness prediction model used is the 
same model that was described in a previous paper. 2 
Previously, this model was not considered reliable 
due to the deteriorating accuracy of the code for mul- 
tiple time steps. Due to the increased accuracy of the 


code for multiple time steps, this routine was reacti- 
vated. Comparisons will be made between this 
model and the measured roughness heights 
obtained by Shin. 3 This routine is considered reliable 
enough that the standard input of sand-grain rough- 
ness into LEWICE has been removed. 

The third feature is the addition of multi-element 
capability to LEWICE. The potential flow solver has 
always been capable of producing multi-element 
flows, but only now have the trajectory, energy bal- 
ance and ice addition routines been correctly modi- 
fied to produce multi-element ice accretions. A 
comparison will be made with experimental data 

obtained in the IRT. 4 

Since several versions of the LEWICE code are 
being used by industry, the next section will provide a 
short history of LEWICE along with the current 
nomenclature of the codes. Results from the modifi- 
cations will then be presented. 

LEWICE History 

In 1983, as a result of university grants and in- 
house research at the NASA Lewis Research Center, 
three computer codes were developed: a potential 
flow code, a droplet trajectory code, and an energy 
balance code. This combined effort, which was called 
LEWICE, was used exclusively for in-house research 
at Lewis. In current nomenclature, this will be called 
LEWICE version 0.1 . 

Through funding by FAA and NASA Lewis, in 
1987 the previous codes were combined into a form 
usable by industry and distribution began. This ver- 
sion will be called LEWICE 0.5. 

Through additional funding by NASA Lewis, 
interactive graphics capabilities were added and a 
correlation for surface roughness were added. The 
code’s capabilities and usefulness were documented 
in CR 185129, Users Manual for the NASA Lewis Ice 
Accretion Code , LEWICE. This version will be called 
LEWICE 1.0 

As usage of the code increased, both in industry 
and at NASA Lewis, several errors were detected 
and fixed in the code and several new features were 
added. These new capabilities were documented in a 
previous AlAA paper 2 and a NASA Contractor 
Report. 5 This version was initially released in June 
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1993 as LEWICE Beta, and will be called LEWICE 
1.3 in current nomenclature. 

The oode described in this paper is the first 
update to LEWICE since that time. It is available for 
limited release. Users of this code, which will be 
called LEWICE 1 .6, are advised that it is still under 
development and may contain several bugs which 
have not yet surfaced. Users who request version 1 .6 
will be asked to promptly report problems they have 
with the code and to provide to NASA the input files 
for these problem cases. 

Future plans for LEWICE upgrades include a 
new roughness model and new heat transfer model, 
both based on experimental data. A new Users Man- 
ual will be produced for that code instead of the cur- 
rent updates. Current nomenclature will refer to this 
future code as LEWICE 2.0 

Numerical Improvements 

Several improvements were made to the 
LEWICE code which allow for more accurate ice 
shape predictions. These include changes to the sur- 
face modeling, mass addition and time stepping 
algorithm to improve convergence characteristics of 
LEWICE. Modifications were also made to the transi- 
tion model, the calculation of transition heat transfer 
transition heat models to produce more realistic heat 
transfer. New additions include an ice growth algo- 
rithm which allows ice growth in arbitrary directions 
and a new “pseudo" surface which features improve 
heat transfer prediction for large glaze shapes. 

Several changes to the LEWICE calculation 
scheme have been made to achieve good conver- 
gence characteristics. A converged solution is one 
that ceases to change with increasing size and or 
increased number of time steps. Any good numerical 
method should have good convergence characteris- 
tics. Convergence is generally controlled by increas- 
ing the number of grid points or time steps within the 
accuracy of the computer. 

For LEWICE, convergence is measured using 
ice shapes and is controlled by the number of time 
steps and grid points. Convergence has occurred 
when an increase in time steps or grid points results 
in no appreciable change in the ice shape. LEWICE 
has had a history of poor convergence characteris- 
tics. The complex, nonlinear nature of ice growth has 
been a severe impediment in producing conver- 


gence. The number of possible time steps have been 
limited to 5 steps for some complex glaze shapes 
and 10 steps for some simple rime shapes. After 
these limits flow code failure and ice shape diver- 
gence were common. These changes have been 
made to LEWICE to overcome the problems result- 
ing in poor convergence and cases have been run 
using as many as 1000 time steps. 

An adaptive grid technique was incorporated into 
LEWICE to allow optimization for each phase of the 
ice accretion process. Surface models are generated 
from a highly refined baseline model using adaptive 
grid techniques. Separate surface models can be 
used for each step in the ice accretion calculation 
depending on the specific accuracy requirements for 
that step. The adaptive grid method is useful in pro- 
ducing surface models with point distributions that 
are sufficiently dense in regions of interest, vary 
smoothly and have a minimum of points, insuring 
computational accuracy and speed. The adaptive 
grid technique involves the use of weighting func- 
tions supplied by the user to produce weighted point 
distributions for a surface. Surface models can be 
generated which resolve any variable of interest in 
the icing process. Weight functions can be chosen to 
produce densely packed points in regions of interest 
and sparsely packed points away from the region of 
interest. Possible weighting functions could include: 
radius of curvature, velocity or velocity gradient, col- 
lection efficiency or collection efficiency gradient, ice 
thickness or ice thickness gradient. 

Currently three surface models with various 
options are used for the icing calculation. For the 
heat transfer calculation two types of models can be 
used: a standard surface model and a “pseudo" sur- 
face model. For both models, point distributions are 
weighted using radius of curvature with user speci- 
fied constraints on maximum and minimum point 
spacing. For the trajectory calculation, surface point 
distributions are optionally weighted to collection effi- 
ciency or towards a constant leading edge spacing. 
The baseline model is used for the remaining mass 
balance, energy balance and the mass addition cal- 
culations. This model, which features several thou- 
sand equally spaced points, was found to be 
sufficiently accurate for these calculations and 
reduced the error in the task of updating the baseline 
model after each time step. 

The ice growth scheme has been modified to 
produce more conservative, accurate, and smoother 
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ice shapes. The scheme involves a separate time 
stepping procedure tor the ice growth module. Ice is 
added to the geometry in small increments to a 
highly refined surface model. During each ice addi- 
tion step, the new surface model is checked for 
regions of convergence and divergence. Points in 
regions of high convergence are removed to keep 
the region from growing into itself. Points are added 
in regions of high divergence to preserve the resolu- 
tion of the surface model. Step size is varied accord- 
ing to convergence and divergence rates of the 
surface model to produce a smooth ice shape. The 
ice addition routine continues, conserving mass 
locally, until all of the ice has been added for the cur- 
rent time step. 

A new variable time stepping algorithm was 
installed which has better convergence characteris- 
tics, is more automated and is computationally more 
accurate and quicker then traditional constant time 
stepping algorithms. The new method involves 
restricting the maximum ice thickness for a given 
time step to a specified amount. Essentially the max- 
imum allowable ice thickness is set and the time step 
is allowed to vary to produce this maximum thick- 
ness. A default value of 0.2% chord has been found 
to be reasonable, but this can be changed by the 
user. 

This type of time stepping allows direct control 
over the most important parameter affecting conver- 
gence which is geometric change. Geometric change 
is the direct cause of change in both collection effi- 
ciency and heat transfer distribution which are the 
two main contributors to ice shape calculation. If the 
amount of geometric change can be controlled in a 
systematic way, then changes in collection efficiency 
and heat transfer distribution can be controlled as 
well and good convergence characteristics can be 
achieved. 

This time stepping technique should be highly 
independent of both geometry and icing condition 
lending itself to automation. Simply put, if the maxi- 
mum geometric change which produces a discernible 
change in collection and heat transfer is known, then 
this value can be set and should not have to be 
changed for any condition or configuration. The tech- 
nique produces larger times steps for periods of 
small geometric growth and smaller time steps for 
periods of large geometric growth than the traditional 
methods, producing fewer and more accurate time 
steps. 


The transition model has been modified to pro- 
duce realistic transition for cases with multipie stag- 
nation points. Previously it was possible for transition 
to be predicted prematurely for cases with multiple 
stagnation points. This resulted in poor heat transfer 
distributions where transition would be predicted very 
near the leading edge. This would cause the heat 
transfer coefficient aft of the transition point to be 
less than the laminar value. This shortcoming is 
responsible for a good portion of the ice shape asym- 
metry observed for seemingly symmetric glaze icing 
conditions (e.g. NACA-0012 at 0 degrees angle-of- 
attack). The transition model is highly dependent on 
velocity gradients, which are realistically large 
between multiple stagnation points. Transition was 
predicted prematurely between consecutive stagna- 
tion points because of these large gradients even 
though the velocities were near zero. The newly 
incorporated model forces transition to be delayed 
until aft of the aftmost stagnation point at the leading 
edge producing more realistic heat transfer. This cor- 
rection provides the desired result of producing tur- 
bulent heat transfer coefficients which are higher 
than the laminar values, which is a more realistic 
result. 

Corrections were made to the transition heat 
transfer calculation so that the boundary layer would 
not prematurely transition from laminar to turbulent. A 
premature transition occurs for high roughness val- 
ues and causes the turbulent heat transfer coefficient 
downstream of transition to be lower than the laminar 
value. Transition is delayed in these cases until the 
computed turbulent heat transfer is greater than or 
equal to the laminar value at this location. 

A new ice growth algorithm which allows specifi- 
cation of arbitrary ice growth directions has been 
implemented. Current options include growth direc- 
tions in the surface normal direction, flow direction 
and in the trajectory tangent directions. This algo- 
rithm will allow easy incorporation of future ice 
growth models. 

Finally, a new optional *pseudo' surface has been 
implemented which produces more realistic heat 
transfer for cases with large stagnation zones than 
the previous model. It is known that inviscid codes do 
not produce realistic surface velocity distributions for 
large concave forward facing regions such as a large 
glaze ice shape. Surface velocity and surface veloc- 
ity gradients are overpredicted resulting in poor heat 
transfer prediction. A more realistic flow solution can 


4 of 12 


be obtained using the inviscid panel codes by filling 
in the large concave region, essentially modelling it 
as a forward facing flat plate. A method has been 
implemented which scans the surface model for 
voids and fills them. 

As evidence of the improved capabilities this 
model provides, the following case was run. A 21 
inch NACA0012 was run at 0° angle of attack for 45 
minutes using 10 second time steps, resulting in 270 
time steps in the simulation. 10 seconds is smaller 
than the smallest time step used in the automated 
time step mode. The automated time step procedure 
used 83 time steps, an average of 32.5 seconds per 
time step. 

Two conditions were run. Both had a velocity of 
150 mph, liquid water content of 0.5 g/m 3 , and a 
droplet diameter of 20 microns. The first case was a 
rime ice case with a total temperature of -8 °F while 
the second case was a glaze ice case with a total 
temperature of 28 °F. 

The final shape of the rime ice case is shown in 
Figure 1 . Although experimental data is not taken for 
this long of an icing time, the mass, shape and maxi- 
mum thickness are proportional to shorter cases ran 
in the IRT to which LEWICE has been compared. To 
demonstrate the high accuracy of the current model, 
the lower surface is reflected upwards to show the 
symmetry of this shape. This is shown in Figure 2. 
The symmetry of the ice is nearly perfect. It was 
unheard of to run the previous versions for more than 
10 time steps, and even then accuracy of this level 
could not be obtained. 

The glaze ice case, ran at a 28 °F total tempera- 
ture, is shown in Figure 3. A 1 0 second time step was 
used for this case as well. Again, although experi- 
mental data is not available for this run, the horn 
growth proceeds logically as time progresses. Previ- 
ous versions of LEWICE could only be run accurately 
for 5 or 6 time steps in glaze conditions, yet this ver- 
sion can easily handle 270 time steps with such large 
horns. Again, the lower surface is reflected upwards 
to show the symmetry of the solution. This result is 
shown in Figure 4. 

The symmetry achieved numerically far exceeds 
the symmetry available experimentally for much 
shorter icing times. The difference between the upper 
and lower horn angles is 0.7° and the difference in 
horn length is 2 mm (4% of length). The symmetry 
achieved here shows that each module is being exe- 


cuted accurately. If any of the modules is not being 
performed correctly, then the shape would not be 
symmetric. This result provides confidence that the 
ice shape is being computed accurately for other 
cases. 

Figure 1 
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Figure 2 

Symmetry Comparison for 45 Minuts les Ships 
270 tims steps of 10 seconds each 



Figure 3 
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Rout# 4 



Roughness Prediction 

The model used for predicting the distribution of 
roughness on an airfoil was developed in a previous 
paper, therefore a detailed derivation will not be pre- 
sented. The height of roughness is determined by 
equating the force of the water flow with the force 
caused by surface tension. The runback model used 
to arrive at these forces was developed by Al-Khalil. 6 

The advantages of this model are that it pro- 
duces results which qualitatively agree with mea- 
sured roughness data and that it relieves the burden 
on the user of having to input a crucial parameter for 
glaze ice. The disadvantage of the current model is 
that it does not predict some of the trends shown by 
the measured roughness data. The equation for the 
local height of roughness is 



where the upper bound on this height is the local 
height of the ice shape. The use of this upper bound 
shows the first discrepancy in the experimental data. 
The measured roughness reached a stable height 
after two minutes and did not increase thereafter, 
whereas the predicted roughness will continue to 
increase. 


An example of the roughness distribution 
obtained by this method is shown in Figure 5. The 
conditions represent the 'baseline' case used by 
Shin 3 in his IRT test. His measured data is repre- 
sented by the squares. The size of the squares rep- 
resents the margin of error in the data. 

Rgur* 5 



In that effort, he reported a ‘smooth’ zone near 
the leading edge, followed by a region of uniform 
roughness. This is qualitatively represented in the 
model, although the predicted distribution does not 
show an abrupt change in height. Shin then reports a 
feather height, which is also shown in this plot 
Despite the agreement of the roughness distribu- 
tion, the distribution itself is not used in the code due 
to the numerical inaccuracy of the flow derivative. 
After several time steps, the quality of the solution 
decreases significantly if the distribution is used. 
Instead, an average roughness is used in the code. 
The equation for this is 


\bids 

*«., = V- (2) 

jds 

o 

This average roughness is then used as the 
sand-grain roughness for the entire airfoil. This aver- 
age value is amazingly close to the measured rough- 
ness values found in Shin’s experiment. In the test, 
parameterizations of LWC, velocity, temperature and 
time were performed to establish the variation of 
roughness with these parameters. The tests were 
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conducted on a 21" chord NACA0012 at 0° angle of 
attack. The matrix centered around these conditions: 

LWC = 0.5 g/m 3 
V = 150 mph 
T = 28 °F 
MVD = 20 pm 
time = 6 min. 

The parameterization used three velocities (150, 
200, 250 mph), three temperatures (25, 28, 30 °F), 
four LWC's (0.5, 0.75, 1.0, 1.2 g/m 3 ) and four times 
(1, 2, 3, 6 minutes). 

The first case shows the increase in roughness 
with icing time and is shown in Figure 6. For this plot, 
the experimental data is shown as solid squares with 
error bars. These error bars are 0.06 mm and repre- 
sent an estimate by Shin of the error in the measur- 
ing technique. He also reported errors for each 
individual measurement, but it was more convenient 
to show a single error level for all measurements. 
The model predicts the level of roughness, although 
the increase with time is more gradual in the pre- 
dicted values. It also continues to increase with time, 
whereas the experimental data shows no significant 
increase in roughness past 2 minutes. However, it is 
an improvement over previous models which had no 
variation in roughness with time. 


Figure 6 

Variation of Roughnasa with Tima 



The variation of roughness with velocity is shown 
in Figure 7. The model shows excellent agreement 
with this parameter. This is an important result, as 
the previous model for sand-grain roughness pre- 
dicted an increase in roughness, as shown by the 
dotted lines. The values are nearly constant over this 
range of conditions due to two competing effects in 
the theoretical model. The increase in velocity 
causes an increase in the amount of incoming water, 


which should increase the roughness. However, the 
higher velocity also causes a change in the flow 
derivative which results in a decrease in the rough- 
ness levels, especially near stagnation. The two 
effects cancel each other, causing essentially no 
change in roughness for these conditions. 


Figure 7 

Variation of Roughnoas will Votodty 



The effect of temperature on roughness is shown 
in Figure 8. The roughness values predicted are rea- 
sonably close to the measured values, however the 
trend is opposite. The measured roughness 
increases very slightly with temperature, while the 
predicted levels show a minimum value at 28 °F, with 
increased values at both ends, especially at the 25 
°F condition. 


Figure 6 

Variation of Rough with Temperature 



The increase in roughness from the 28 °F case 
to the 30 °F case is caused by a slight increase in 
surface tension with temperature. The increase in 
roughness from the 28°F case to the 25 °F case is 
caused by a decrease in the “wetting factor’ which is 
caused by a change in the predicted contact angle of 
the roughness element. This theory needs to be 
developed further to correct the predicted trend. The 
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current roughness prediction and the previous corre- 
lation are both within the experimental error over this 
range. 

The variation of roughness with liquid water con- 
tent (LWC) is shown in Figure 9. The roughness 
increases with LWC in both the predicted and mea- 
sured values, but there is a significant (20%) under- 
prediction in roughness at the higher LWC values. 
However, this model is a vast improvement over pre- 
vious models which predicted extremely high rough- 
ness values at high LWC, as shown by the dotted 
lines. 


Figure 0 

Variation of Roughnaas with LWC 



Multi-Element Ice Accretion 

The Hess-Smith potential flow code used by 
LEWICE had the capability to predict flow on multi- 
element configurations. Often, the limitations of 
potential flow make these predictions much less 
accurate than more sophisticated models, but it may 
be adequate for the purpose of providing sample ice 
shapes. However, the trajectory and ice accretion 
routines were only applicable for single element 
geometries. As a result of the current effort, not only 
can LEWICE perform multi-element trajectories, but 
these routines can also be converted for use with 
more accurate flow solvers. 

LEWICE performs multi-element trajectories by 
treating each element as a single entity. Impinge- 
ment limits and collection efficiencies are determined 
on each body as though the other bodies are not 
there. Their influence on the trajectories is embed- 
ded in the flow solution, which takes into account all 
of the bodies. Any trajectory hits on other elements 
are treated as missed trajectories. 


The hits on other bodies are, however, useful in 
determining the starting location of the next trajectory 
in the impingement limit search, especially hits on 
bodies which precede the one selected. For exam- 
ple, when the code looks for impingement limits on 
the flap(s), trajectory hits on the slat are useful in 
determining the starting location of the next trajec- 
tory. Routine MODE in LEWICE determines if a tra- 
jectory hits or misses a body. It was modified so that 
it not only knew that a body was hit, but which one. 
As stated earlier, hits on other bodies are only used 
to select the starting location of the next trajectory. 

An additional problem occurs e specialty on the 
main element of a multi-element airfoil. It is possible 
for trajectories to hit this element by passing above 
the slat as well as by passing below the slat. There- 
fore, for all bodies but the first one (the slat) LEWICE 
will first look below the slat for an upper and lower 
impingement limit and determine one set of collection 
efficiencies for this set of impingement limits. 
LEWICE will then look above the slat and attempt to 
find a second set of impingement limits. If two sets of 
limits are found, the two collection efficiency arrays 
are merged. 

This process is made clearer by looking at Fig- 
ures 10 and 11. Figure 10 shows the two sets of 
impingement limits on the main element of a sample 
slat and main element combination. The collection 
efficiency curve for this condition is shown in Figure 
11. Both bodies are NACA0012 airfoils. This is not a 
realistic test case, but is simply representative of the 
code’s capabilities. If the user knows that impinge- 
ment is going to occur on other elements by travel- 
ling below the slat, the user can bypass this option as 
the code will run nearly twice as fast by bypassing 
this feature. 


Figure 10 
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Figure 11 

Collection efficiency for 2 element test ceee 



The process of converting LEWICE to handle 
multi-element geometries was made more difficult by 
the addition of the multiple geometry scheme used to 
increase the accuracy of the code. After the trajectory 
routine is completed, LEWICE creates a second set 
of panels for each element and resolves the flow field 
using this panel set. For a single element case, this 
second set of panels produces a smoother pressure 
distribution which increases the accuracy of the 
boundary layer integration. This is not always the 
case for multiple element geometries. Work is con- 
tinuing on this routine so that this very useful feature 
can be used with multi-element geometries. 

Once the collection efficiencies for each element 
has been found, and the flow recalculated, LEWICE 
is ready to perform the boundary layer integration 
and the ice accretion. Once again, this procedure is 
performed on each element individually, without 
regard to the presence or influence of other ele- 
ments. The geometry modification is performed on 
each element individually, hence the code does not 
currently check for different elements intersecting 
due to ice growth. In this case, the code will most 
likely crash when it tries to solve the flow field on the 
next time step. 

A comparison between this code and experimen- 
tal data taken on a 5-element Boeing 737-200 airfoil. 
The airfoil is shown in Figure 12. The experimental 
data was taken in the IRT in 1 991 and is documented 
in reference (4). The conditions for the comparison 
were: 

15° flap 
V = 100 mph 
T = 28 °F 


LWC = 0.92 g/m 3 
MVD = 14.4 pm 
0° angle of attack 
time = 8 minutes 


Figure 12 



This case was selected because it was a glaze 
ice condition and ice was obtained on all five ele- 
ments. The ice shape comparison is shown in Figs. 
13-17. Due to the complex geometry, a Langmuir ‘D’ 
droplet distribution consisting of 7 drop sizes was 
used in the numerical prediction. 


Figure 13 
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Figure 14 
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Figure 15 
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Figure 17 



extensive runback. This prediction is quite far from 
the experimental icing limits. 

There are several reasons why this occurred. 
First, there are questions as to whether the configu- 
ration used in LEWICE is the same as the test config- 
uration. Airfoil coordinates are only available in the 
stowed reference plane. A rotation angle, x-offset, 
and gap distance are listed in a table in the report 
which allows coordinates for the test configuration to 
be computed. However, the 15° flap configuration is 
the only one of the three configurations for which the 
elements do not intersect with each other when this 
transformation is applied. This leaves open the ques- 
tion as to the validity of all of the numbers in the 
table, including this configuration. By observing pic- 
tures of the airfoil in the 15° flap configuration, it is 
felt that the correct configuration was tested with 
LEWICE. 

The second explanation for this discrepancy is 
the quafity of the potential flow solution, which is plot- 
ted in Figure 18. The flow solution may not be ade- 
quate for this complex flow situation. A resolution of 
this question will not be known until more derailed 
analysis can be made with more accurate flow solv- 
ers. 


Figure 15 
Predicted Pretture Coefficient 



The predicted ice shape on the slat, shown in 
Figure 13 is representative of the very poor predic- 
tion when compared to the experimental data. This 
poor performance is somewhat surprising, as previ- 
ous improvements to LEWICE had shown significant 
improvement in ice shape prediction for single ele- 
ment airfoils. The extent of predicted ice accretion in 
this figure is due to direct impingement, not due to 
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Figure 19 

Particle T r ajec t ories tor a 



Finally, multiple time steps of this case were not 
possible due to errors within the code. Currently, only 
single time step ice shapes are possible for multi-ele- 
ment configurations until this issue is resolved. The 
changing iced geometry can change the collection 
efficiency distribution on a single element airfoil dras- 
tically and will no doubt have an impact on the results 
for this case. 

The other elements also show major discrepan- 
cies between the predicted ice shape and the experi- 
mental ice shape. The ice shape for the main 
element is shown in Figure 14. In both the numerical 
and experimental shapes, a small amount of ice is 
produced just past the slat. The predicted shape is 
further down the chord than the experimental shape. 
This is consistent with the difference shown in the 
slat ice shapes. If the slat is raised slightly in space, 
more water will hit the main element and the 
impingement will be closer to the leading edge. It is 
unknown if this is due to errors in the flow or an error 
in the definition of the slat geometry. 

The comparison for the fore flap is shown in Fig- 
ure 15. The predicted impingement limits for this ele- 
ment are also much wider than the experimental 
shape. The prediction for this element would also be 
greatly improved if additional time steps could be ran. 
The prediction for the other two elements show the 
same trends. The prediction for the main flap, shown 
in Figure 16, shows some significant runback freez- 
ing. The upper surface shape is due to runback 
freezing and intersects the fore flap, meaning that the 
gap between them has filled with ice. This effect is 
one of the problems which must be dealt with for 
multi-time step runs to be run. The comparison for 
the aft flap shows the same trends as described ear- 


lier and are likely due to accuracy of the flow solution 
and use of a single icing time step. 

Conclusions 

Three major improvements have been made to 
the LEW1CE code. The LEWICE code has been 
modified to handle a large number of time steps for 
single element geometries and to perform the ice 
addition and recalculation of the flow accurately. A 
roughness model has been proposed which predicts 
measured roughness values much more accurately 
than previous models. Finally, multi-element trajecto- 
ries and single time-step ice accretions on multi-ele- 
ment geometries are now possible. 

Due to the ability of the code to handle a great 
number of time steps very accurately, little future 
work will be performed in this area due to the suc- 
cess of this routine. Work on the multi-element fea- 
tures is needed so that multiple time steps can be 
performed on these configurations as well. 

The roughness prediction model shows much 
promise, but additional work is needed on a more 
robust model. Research is also being needed on the 
affect of this roughness on the heat transfer. 

The code has shown the ability to predict collec- 
tion efficiencies and single time-step ice accretions 
on multi-element geometries. Additional work is 
needed to extend this routine to multiple time steps 
and to make the solution more robust. 

This new version of LEWICE is called LEWICE 
1 .6 (Beta) and is available for use by industry. Indus- 
try participants are forewarned that development of 
this code is continuing, especially for the multi-ele- 
ment features. For single element geometries, the 
code is working well. The multi-element option is still 
under development and may contain ‘bugs’ which 
have not yet been discovered. Care must be taken 
when using this option. Limited release of version 1 .6 
is recommended for multi-element users until these 
problems are resolved. The multi-element capability 
for LEWICE is considered a transitional code until 
these multi-element trajectory capabilities can be 
added to better multi-element flow codes. 
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